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ABSTRACT 

It has already been shown, using a local model, that 
accretion discs with cooling times icooi < 351"^ frag- 
ment into gravitationally bound objects, while those 
with cooling times icooi > 351^^ evolve into a quasi- 
steady state. We present results of three-dimensional 
simulations that test if the local result still holds 
globally. We find that the fragmentation boundary 
is close to that determined using the local model, 
but that fragmentation may occur for longer cool- 
ing times when the disc is more massive or when 
the mass is distributed in such a way as to make 
a particular region of the disc more susceptible to 
the growth of the gravitational instability. These re- 
sults have significant implications for the formation 
of of gaseous planets in protoplanetary discs and also 
for the redistribution of angular momentum which 
could be driven by the presence of relatively mas- 
sive, bound objects within the disc. 



1. INTRODUCTION 



The discovery of th e first extra-solar planet 
l|Mavor & Ouelod . 11995') and the subsequent discov- 
ery of many (> 100) additional extra-solar planets 
(iMarcv fc Butleil 2000) has enhanced the interest in 
the formation of planets and the evolution of proto- 
planetary d iscs. The observation of the first transit- 
ingplanet ijHenrv et all l200nt ICharbonneau et all 
I2OOOD . giving an estimate of the planet radius, and 
the relatively large masses (> O.lMjupitcr) of all 
the currently detected extra-solar planets, has led 
to the view that they are all giant gaseous plan- 
ets. The most widely studied, and accepted, for- 
mation mechanism for gaseous planets is core ac- 
cretion ijLissaueil Il993j^ in which planetesimals grow 
by direct collision to form a core which, when siffi- 
ciently massive (m ^ lOmoarth), then accretes an en- 
velope of gas from the disc. Models l)Pollack et all 
^996), however, suggest that the timescale for planet 
formation via this mechanism may be longer than 
the lifetime of most protoplanetary accretion disc 



ijHaisch. Lada fc Ladal l2001|) . This has led to a re- 
newed interest in the possibility that gas giant plan- 
ets may have formed directly via a disc ins t abilit y 
iKuineTllT9RTllB^IT99ll200?il:lM^^ et a,lll2f)nl . 

A protoplanetary disc that is sufficiently massive and 
cool may form gravitationally bound gaseous ob- 
jects via the gravitational instability. This mech- 
anism is extremely rapid and differs from core ac- 
cretion in that a rocky core is not required. The 
possibility that this mechanism may play a role 
in gaseous planet formation has been enhanced by 
recent models of the interiors of gaseous planets 
ijGuillot. Gautier fc H ubbard. 1997; GuiUot, 1999) 
which suggest that Jupiter could have a relatively 
small core with mcorc < lOrTioarth- A Keplerian ac- 
cretion disc with sound speed c^, surface density S, 
and epicyclic fr equency k will become gravitationally 
unstable if the lToomrd ©64} Q parameter 



is of order unity. A gravitationally unstable disc 
can either fragment into one or more gravitation- 
ally bound objects, or it can evolve into a quasi- 
steady state in which gravitational instabilities lead 
to the outward transport of angular momentum. 
The exact outcome depends on the rate at which 
the disc heats up (through the dissipation of tur- 
bulence and gravitational instabilities) and the rate 
at which the disc cools. It h as been suggested 
(jGoldreich fc Lvnden-Bell Il96,'^^ that a feedback 
loop may exist. When Q is large the disc is sta- 
ble and cooling dominates, driving the disc towards 
instability. When Q becomes sufficiently small, heat- 
ing through viscous dissipation dominates, and the 
disc is returned to a state of marginal stability. In 
this way Q is maintained as a value of ^ 1. 

It has, however , been shown using a local model 
l|GammieL l2001l) that a quasi-stable state can only 
be maintained if the cooling time icooi > 3r2~^ where 
O is the local angular frequency. For shorter cooling 
times the disc fragments. This is consistent with 
iPickett et all l|1998l. 1200(11) that 'almost isothermal' 
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conditions are necessary for fragmentation, and de- 
fines a robust lower limit to the critical cooling time 
below which fragmentation occurs. It has been sug- 
gested, however, that self-gravitating discs req uire 
a str ictly global treatment ( Balbus & Papaloi zoii 
II999D . and while global effects are highly unlikely 
to stabilize a locally unstable disc, they could well 
allow fragmentation within discs that would be lo- 
cally stable. This has led to an interest in global 
simu lations of self-gravitating protoplanetary discs 
(Ric e et al.L |^Q3J) and we discuss, in this paper, re- 
sults of a number of these simulations. 



2. NUMERICAL SIMULATIONS 

2.1. Smoothed particle hydrodynamics 

The three dimensional, global simulations presented 
here were performed using smoothed particle hy- 
drodynamics (^PH)j_a^Jjagrangiaji hydrodynamics 
code (|Benzl 1199(1 iMonaghani 119921) . The central 
star is modelled as a point mass that may accrete 
gas particles if they approach within a predefined 
sink radius, while the gaseous disc is modelled us- 
ing 250000 SPH particles. A tree is used to de- 
termine neighbours and to calculate gravitational 
forces between gas particles and between gas par- 
ticles and point masses (j3cnz, 199(|). If in the sim- 
ulation the disc starts to fragment, producing high 
density regions, the code can in principle continue. 
The high density regions do, however, slow the code 
down significantly. To continue simulating a frag- 
menting disc, high density regions that are grav- 
itationall y bound are c onve rted into point masses 
( Bate. B onnell fc PricS [l99.5) which may continue 
to accrete disc gas. Since the likely number of point 
masses is quite small, compared to the number of 
gas particles, the gravitational force between point 
masses is computed directly. An additional saving 
in computational time is also made by using individ- 
ual particle time-s t eps JBate. Bonnell fc Pric3 . ll99fit 
iNavarro fc Whitd . Il993j) . The time-steps for each 
particle is limited by the Courant condition and by 
a force condition ijMonaghanl Il992|) . 



2.2. Initial conditions 

We consider a system comprising a central star, mod- 
elled as a point mass with mass Af*, surrounded by 
a gaseous circumstellar disc with mass Afdisc- We 
performed a number of simulations with disc masses 
of Mdisc = O.IM* and Mdisc = 0.25M*. The disc 
temperature i s taken to have an initia l radia l profile 
of T cx r-°-^ l|Yorke & Bodenheimerl [1399?) and in 
most of the simulations the Toomre Q parameter is 
assumed to be initially constant and to have a value 
of 2. A stable accretion disc where self-gravity leads 
to a steady outward transport of angular momen- 
tum should have a near constant Q throughout. A 
constant Q together with Equation 1 gives a surface 



density profile of E oc r~^/^, and hydrostatic equilib- 
rium then gives a central density profile of p oc r~^. 
We did, however, perform a single simulation with 
Afdisc = O.IM* and with S oc . In this case, 
T c>c r~^^^ gives a Q that is not constant but that 
decreases with increasing radial distance. The tem- 
perature was chosen such that the initial minimum 
value of Q was 1.5. 

The disc is modelled using 250000 SPH particles, 
which are initially randomly distributed such as to 
give the specified density profile between inner and 
outer radii of rin and rout- Since the calculations are 
scale free, we take = 1, nn = 1, and rout = 25. If 
we were to assume a physical mass scale of 1 Mq and 
a length scale of 1 au, the central star would have a 
mass of 1 Mq, the circumstellar disc would have a 
mass of 0.1 Mq or 0.25 Mq and would extend from 
1 to 25 au, and 1 yr would equal 2n code units. 



2.3. Cooling 

The main a im of thi s work is to test if the results 
obtained bv lGammi i pOOH) still hold globally. Con- 
sequently we use the same approach in these sim- 
ulations as was used in the local model. We use 
an adiabatic equation of state, with adiabatic index 
7 = 5/3, and allow the gas to heat due to both PdV 
work and dissipation. Cooling is implemented by 
adding a simple cooling term to the energy equation. 
Specifically, for a particle with internal energy per 
unit mass Ui, 

dt ^cool 

where, as in ICamrni^ (|2001l) . tcooi is given by /3rt~^ 
with the value of (3 varied for each run. 

Although the above cooling time is essentially chosen 
to compare the local model results with results using 
a global model, it can also be related (at least ap- 
proximately) to the real physics of an accretion disc. 
For an optically thick disc in equilibrium, the cool- 
ing time is given by the ratio of the thermal energy 
per unit area t o the radiative losses per unit area. It 
can be shown (jPringld I198H) that in such a viscous 
accretion disc, the cooling time is 

97(7-1)^ 

wh ere fl is the angular frequenc y in the disc, and a is 
the lShakura fc SunvaevI (^2^ viscosity parameter. 



3. RESULTS 



3.1. Mdi.c = 0.1M, 

We use a three dimensional, global model to consider 
how cooling times of icooi = 50~^, and ^cooi = 
affect the gravitational stability of a protoplanetary 



Figure 1. Equatorial density structure for a disc with 
-^disc = O.IM*, S oc r"^/* and with a cooling time of 
icooi = 5f2^^. The density is highly structured with 
the instability present at all radii. The density has, 
however, not increased significantly, and the disc is 
in a quasi-steady state with the heating through the 
dissipation of the instability balancing the imposed 
cooling. 



accretion disc with Mdisc — O.IM.^ and with S oc 
j,-7/4^ This choice of surface density profile, together 
with T (X r""'^, gives a Toomre Q parameter that is 
intially constant throughout the disc. 

Figure n] shows the final equatorial density structure 
of the icooi = 5f2^^ simulation. The central star (not 
shown) is located in the middle of the figure, and 
the X and y axes both run from —25 to 25. The 
disc is highly structured and the instability exists at 
all radii. However, at no point in the disc has the 
density increased significantly, and no fragmentation 
has taken place. Figure El shows the Toomre Q pa- 
rameter for the same simulation at three different 
times during the simulation. At the beginning of the 
simulation (t = 0) Q has an almost constant value 
of 2. At the end of the simulation {t = 2932) the 
value of Q is almost unity between radii of 1 and 15. 
Comparing Q at i = 876 and t = 2932 shows that 
the cooling initially reduces Q to a value below 1, 
causing the instability to grow, heating the disc and 
returning Q to a value of order unity. The disc has 
clearly settled into a quasi-steady state in which the 
imposed cooling is balanced by heating through the 
dissipation of the gravitational instability. 

Figure shows the equatorial disc structure for a 
simulation with the same disc mass and surface den- 
sity profile as in Figure^ but with a cooling time of 
icooi = 30""'^. This simulation was only run for 504 
time units and hence we show only the inner 8 radii 
of the simulation. The disc is again highly struc- 
tured, but in this case the disc has started to frag- 
ment into gravitationally bound clumps, indicated 
by the bright dots in Figure |3 The density in the 
clumps is 4 - 5 orders of magnitude greater than the 
initial disc density, and to reach the time shown in 



Figure 2. Toomre Q parameter at the beginning (t = 
), one third of the way (t — 876 ), and at the end 
(t = 2932 ) of the simulation in which Mdisc = O.lAf,, 
S oc r~'^l^, and ^cooi = 551""'^. At the end of the 
simulation Q is of order unity between radii of 1 and 
15. 

Figure 01 many of the clumps have beeri converted 
into point masses ijBate. Bonnell fc Pried Il995j) . 

Figure 01 shows the Toomre Q parameter at three dif- 
ferent times during the above simulation. At all three 
times shown, there is a region which has Q signifi- 
cantly less than 1. In the regions where Q is less than 
unity, the instability grows rapidly and produces 
gravitationally bound fragments. Once formed, these 
fragments tidally heat the disc, increasing Q to a 
value greater than unity. Once Q > 1, fragmentation 
ceases and the disc, at that radius, becomes gravi- 
tationally stable. There is, however, likely to be a 
significant amount of angular momentum transport 
driven by tidal interactions between the gravitation- 
ally bound fragments and the gaseous disc (,Larsork 
l200l . 



3.2. Mdisc = 0.25M* 

To study how disc mass may influence the global na- 
ture of the gravitational instability we consider a disc 
with a mass of Mdisc = 0.25M*. As in the previous 
simulations, the surface density is taken to have a 
radial profile of S oc which, together with T oc 

^-o.b ^ gives an intially constant Toomre Q parame- 
ter. Although most T Tauri discs have masses con - 
siderably less than 0.25Mh. l|Beckwith et all 1199(H) . 
there a few with such masses, and this simulation 
may also apply to an earlier stage of the star for- 
mation process when discs are expected to be more 
massive. 

Figure El shows the final equatorial density structure 
of the above simulation. The disc is highly struc- 
tured and the nature of the spirals, compared to the 
equivalent simulation with Mdisc = .1M<,, is con- 
sistent with the increased disc mass llNelson et al.L 



Figure 5. Equatorial density structure of a disc with 
Mdisc = 0.25M,, E oc r-y^, and tcooi = ■ 
There are clear signs of fragmentation with the most 
massive fragments being gravitationally boud. 



Figure 3. Equatorial density structure of the inner 
regions of a disc with Mdisc = O.lAf,, E oc r"^/**, 
and with a cooling time of tcooi = 3f2~^. The disc 
is highly structured and is starting to fragment into 
high density, gravitationally bound clumps. 




Figure 4- Toomre Q parameter at times of t = 192, 
t = 392, and t = 504 for Mdisc = O.IM*, S oc r^^/^, 
and icooi = 3f^^^- 



Il998|) . Unlike the equivalent lower mass simulation 
(see Figure ^| , there are a number of high density 
clumps present in the disc. The standard routine for 
checking if these clumps are bound initially found 
them to be unbound. This routine, however, only 
considers the nearest ~ 50 SPH neighbours. By in- 
creasing this to the nearest ~ 250 SPH neighbours, 
the densest clump was found to be just gravitation- 
ally bound. This would suggest that in more mas- 
sive discs, global effects could act to make the disc 
more unstable, allowing gravitationally bound frag- 
ments to grow for cooling times greater than that 
obtained using a local model. This could have signif- 
icant implications for the formation of planets, or bi- 
nary companions (Adams, Rudcn & Shu, 1989), rea- 
sonably early in the star formation process. 



3.3. 



O.IM*, E oc r" 



All of the previous simulations were performed as- 
suming E oc r^^/**, and T oc r~°'^, giving an ini- 
tially constant Toomre Q parameter. We have per- 
formed a single simulation with A^disc — O.IM*, 
T oc r~'^^^ , and E oc r~^ . With these parameters, the 
Q value decreases with increasing radius. We there- 
fore normalised our temperature such that Q = 1.5 
at r = 25. Figure El shows the final equatorial den- 
sity structure of a simulation with the above disc 
parameters, and with an imposed cooling time of 
^cooi = 5ri~^. The disc is again highly structured 
and unlike the simulation with the same disc mass 
and cooling time, but with the steeper surface den- 
sity profile, there is clear evidence of fragmentation 
in the outer regions of the disc. The shallower sur- 
face density profile (E oc r~^) means that, compared 
to the steeper surface density profile (E oc r^^/**), 
there is more mass at large radii. This suggests that 
the global nature of the instability may depend both 
on the mass of the disc, and on how the mass is dis- 
tributed within the disc. 
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Figure 6. Equatorial density structure of a disc 
with Mdisc = O.lil/,, E oc r~^ , and tcooi — Sfi"^^. 
There are clear signs of fragmentation suggesting 
that the fragmentation boundary may depend both on 
the mass of the disc and on how the mass is dis- 
tributed. 

4. CONCLUSIONS 

It has been shown using a local model 
I2OOI) that a disc will fragment for cooling times 
icooi < Sri"^ and will settle into a quasi-steady state 
for coohng times tcooi > Sf^"^. We present here re- 
sults from three-dimensional simulations which test 
if the local results still holds globally ijRice et all 
12003'). We inrpose the same cooling function as in 
ICammic (200^ which, alth ough fairly sim plistic, can 
also be justified physically l|Pringlel llQSlI^ . 

For a disc mass of Mdisc = O.IM*, a coohng time 
of fcooi = 5r2~^, and an initially constant Q, the 
disc settles into a quasi-steady state in which Q is 
of order unity over a large region of the disc. De- 
creasing the cooling time to icooi = 30~^ causes the 
disc to rapidly become unstable and produces numer- 
ous gravitationally bound clumps. For this disc mass 
and surface density profile (S oc r^^/'*), the fragmen- 
tation boundary therefore seems to agree with that 
obtained using a local model. 

Increasing the disc mass to Mdisc = 0.25M* and 
keeping all other parameters unchanged, however, re- 
sults in fragmentation for a cooling time of tcooi — 
A similar result was obtained for a simulation 
in which the disc mass was kept at A/disc = O.IM* 
but in which the surface density profile was changed 
from E cx r^''/'* to E oc r~^ . These results suggest 
that global effects may become significant for mas- 
sive disc or for discs in which the mass is distributed 
in such a way as to make a particular region of the 
disc susceptible to the growth of the gravitational 
instability. For T Tauri discs, which generally have 
masses Mdisc < O.IM* ijBeckwith et al.l Il990() . the 
fragmentation boundary is therefore likely to be close 
to that determined using the local model. Earlier in 
the star formation process, when disc may be more 



massive, fragmentation may occur for longer cooling 
times than that predicted by the local model results 
iCammieLim. 

For qu iescent T Tauri discs, the lShakura fc SunvaevI 
l)1973f) a is conventially e stimated to be of the 
order of 10'^ or smaller l|Hartmann et all Il998t 
iBell fc L in. 1994). The simple estimate quoted ear- 
lier (see Eq. |2l would suggest that such discs are 
comfortably stable. Boss (2002), however, has shown 
using an approximate treatment of disc heating and 
cooling, that there may be periods when the cooling 
time is comparable to the orbital period. Our results 
suggest that if the disc is fairly massive, such short 
cooling times may open a window of opportunity 
for the formation of substellar objects, probably in 
the fo rm of a multiple system l|Armitage &: Hansenl 
Il999fl . This has important implifications for the for- 
mation of giant gaseous planets in protoplanetary 
discs and may also pla y a role in the rapid transfer 
of angular momentum l)Larsonl |2002|) . 
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